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RXR signaling is predicted to have a major impact in macrophages, but neither the biological consequence nor the 
genomic basis of its ligand activation is known. Comprehensive genome-wide studies were carried out to map 
liganded RXR-mediated transcriptional changes, active binding sites, and cistromic interactions in the context of 
the macrophage genome architecture. The macrophage RXR cistrome has 5200 genomic binding sites, which are 
not impacted by ligand. Active enhancers are characterized by PU.l binding, an increase of enhancer RNA, and 
P300 recruitment. Using these features, 387 liganded RXR-bound enhancers were linked to 226 genes, which 
predominantly reside in CTCF/cohesin-limited functional domains. These findings were molecularly validated 
using chromosome conformation capture (3C) and 3C combined with sequencing (3C-seq), and we show that 
selected long-range enhancers communicate with promoters via stable or RXR-induced loops and that some of the 
enhancers interact with each other, forming an interchromosomal network. A set of angiogenic genes, including 
Vegfa, has liganded RXR-controlled enhancers and provides the macrophage with a novel inducible program. 
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RXR is an unusual and somewhat neglected member of 
the nuclear receptor superfamily because it is still not 
known how this receptor interacts with the genome and 
regulates gene expression upon ligand activation. It is 
expressed in every cell type and is required for postnatal 
life in mice (Mangelsdorf et al. 1992; for review, see Szanto 
et al. 2004). Its presumed main molecular function is to 
regulate the activity of a dozen or so nuclear receptors. 
There is also evidence that it can form homodimers and/or 
have heterodimer-independent signaling capacity (Zhang 
et al. 1992; Szeles et al. 2010). A key concept regarding RXR 
signaling is the permissiveness/nonpermissiveness mutu- 
ally exclusive dual paradigm. According to this, in certain 
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heterodimers such as RXR:PPAR and RXR:LXR, ligand 
activation of RXR results in transcriptional activation,- 
hence, these are permissive heterodimers, while in other 
heterodimers such as RAR:RXR, TR:RXR, and VDR:RXR, 
RXR is suppressed or " subordinated/ 7 and therefore these 
so-called nonpermissive heterodimers cannot be activated 
from the RXR side (Germain et al. 2002). Therefore, the 
activation of all permissive heterodimers present in a par- 
ticular cell type might lead to pleiotropic gene activation 
and engagement of potentially conflicting pathways. 

The existence of pleiotropy, the role of RXR activation 
of permissive heterodimers, and the presence and activity 
of RXR homodimers have been debated and remain 
largely unresolved. 

The fact that certain natural lipids — such as 9-cis 
retinoic acid (RA), docosahexanoic acid, and phytanic 
acid — are able to activate RXR gives further support to 
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the biological role of RXR activation in vivo (for a review, 
see Szanto et al. 2004). There are also potent and selective 
synthetic compounds such as bexarotone (LG 10069) and 
LG100268 (LG268) (Boehm et al. 1994, 1995) that have 
been used to dissect the role of the receptor in various 
biological systems and/or used in therapies. In macro- 
phages, there are several heterodimeric receptors with key 
cellular roles, such as PPAR7 regulating oxLDL uptake and 
processing, LXR regulating cholesterol efflux and immune 
function, and NR4A1 (NUR77) regulating inflammatory 
response. These heterodimeric receptors have been linked 
to the development of atherosclerosis and also immune 
function and provide means to reprogram macrophages (for 
reviews, see Calkin and Tontonoz 2012; Nagy et al. 2012). 
Therefore, it is important to understand how activation of 
RXR contributes to these pathways and potentially to novel 
ones and regulates gene expression in this cell type. 

A key issue in understanding any signal-specific tran- 
scription factor is the determination of the genomic sites 
to which it binds and to link those to the target tran- 
scripts. Recent advances in genome-wide approaches 
allow global assessment of histone modifications and 
transcription factor genomic binding sites (Ostuni et al. 
2013). It is done using microarrays or RNA sequencing 
(RNA-seq) (Mortazavi et al. 2008). However, the temporal 
changes in steady-state mRNA levels diverge from changes 
in transcription,- hence, simply determining steady-state 
mRNA levels will not allow easy identification of primary 
transcriptional events and will not separate those from 
secondary and tertiary ones (Bhatt et al. 2012). The de- 
velopment of genome-wide localization studies, primarily 
ChlP-seq (chromatin immunoprecipitation [ChIP] com- 
bined with deep sequencing) approaches, aided the de- 
termination of transcription factor-binding sites within 
the context of the chromatin structure on the genomic 
scale. Such genome-wide studies suggest that enhancers 
and promoters exhibit distinct chromatin "signatures." 
The characteristic signature for enhancers consists of 
monomethylation of histone H3 Lys4 (H3K4mel), acetyla- 
tion of histones (H3K9ac and H3K27ac), and binding of the 
acetyltransf erase P300 (Heintzman et al. 2007; Koch et al. 
2007; Visel et al. 2009). 

In spite of all these developments, the reliable linking 
of activated enhancers to the regulated gene is still not 
solved, partly due to the lack of markers linking the ge- 
nomic binding site to the regulated transcript and partly 
due to the potentially large distances between the en- 
hancer and the regulated gene. Moreover, the identifi- 
cation of enhancers is even more cumbersome, if not 
impossible, in the case of signal-dependent transcription 
factors (such as RXR), which may constitutively occupy 
genomic sites. In our studies, we tried to solve these 
issues by combining RNA-seq, ChlP-seq, GRO-seq (global 
run-on sequencing), and 3C-seq (chromosome conforma- 
tion capture [3C] combined with sequencing) in a highly 
integrated way to unravel the mechanism of RXR-induced 
transcriptional events in mouse bone marrow-derived 
macrophages (BMDMs) and, as the result of the process, 
discovered and validated a novel biological activity pro- 
moted by the receptor. 



Results 

The transcriptional consequences of RXR activation 
in murine macrophages 

We set out to systematically determine the genomic 
events following RXR ligand activation in BMDMs (Fig. 
1A). The synthetic and selective RXR ligand LG268 was 
applied throughout the studies in a 100 nM concentration. 
We determined the changing transcripts by RNA-seq 
(analyzed by the pipeline shown in Supplemental Fig. 
SI A) and found that selective activation of RXR affects 
the steady-state mRNA levels of hundreds of genes, as 
expected. A hierarchical clustering of the top 200 changing 
genes as a time course is shown in Figure IB. In order to 
identify primary RXR targets, we embarked on analyses to 
determine both the genomic binding events of RXR (RXR 
cistrome) and the enhancers responsible for the transcrip- 
tional changes. 

Determination of the RXR, PU.l, and P300 cistromes 
and mapping active histone marks in BMDMs, 
the effect of ligand 

By carrying out ChlP-seq experiments, we determined 
the cistromes of RXR in the absence and presence of its 
ligand, LG268, and its relationship to the lineage-specific 
transcription factor PU.l and the cofactor P300 along 
with a marker of transcription initiation, H3K4me3, and 
active histone marks, H3K27ac, H4ac, and H3K4me2, 
using the algorithm depicted in Supplemental Figure SIB. 
We determined —5200 RXR genomic binding regions. 
Importantly, this cistrome is not impacted greatly by 
ligand treatment in 60 min. PU.l occupied the highest 
number of peaks (—30,000), with only a minimal rear- 
rangement upon ligand exposure (Supplemental Fig. SIC). 
Although the number of RXR-binding regions does not 
change much upon ligand exposure, peaks gain —30% more 
reads, suggesting that RXR enrichment on the genomic 
regions is enhanced (Fig. 1C). We confirmed these observa- 
tions by ChIP and real-time quantitative PCR (RT-qPCR) on 
binding regions of known directly regulated genes (Fig. ID; 
Nagy et al. 1996; Kersten et al. 2000; Repa et al. 2000). 
Intriguingly, P300 binding follows the genome-wide RXR 
enrichment, suggesting that P300 is likely to be recruited to 
the genome upon RXR activation (Supplemental Fig. SIC). 

The motif rank order under the detected RXR peaks 
reported that PU. 1 had the most enriched motif, followed 
by various combinations of repeats, including DR1 and 
DR4, of the nuclear receptor-binding (half) site. In addi- 
tion, with somewhat lower abundance, two other mac- 
rophage-associated motifs were also detected: C/EBP and 
AP-1 (Fig. IE). As far as the genomic distribution of the 
detected peaks is concerned, —90% were found outside of 
the potentially directly regulated 823 genes identified by 
RNA-seq experiments (i.e., not within 10 kb of their 
transcription start sites [TSSs]) (Fig. IF), suggesting that 
linking binding sites to regulated genes based simply on 
proximity will be difficult. 

Next, we continued to interrogate the cistromes of 
PU.l, RXR, and P300 in the ligand-activated state. The 
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Figure 1. Mapping the transcriptional consequences of selective activation of RXR in murine macrophages. (A) The scheme of 
macrophage differentiation and the chemical structure of LG268. [B) Heat map representation of the expression of genes regulated by 
LG268. Differentiated macrophages were treated with 100 nM LG268 for the indicated time, and total RNA was subjected to RNA-seq 
gene expression analysis (Supplemental Fig. SI A). Biological replicates had at least 70 million reads each. The presented gene set was 
filtered at FDR < 0.1 and a more than twofold change. (C) Average read distribution of RXR peaks in the presence or absence of LG268. 
ChlP-seq data were analyzed as shown in Supplemental Figure SIB (peaks were predicted by MACS2, and the consensus of the 2X2 
replicates was determined by DiffBind); the histogram was generated by HOMER. (#) Normalized throughout the study. (D) RXR 
binding on the indicated individual enhancers confirmed using ChlP-RT-qPCR. Macrophages were treated with 100 nM LG268 for 1 h. 
The mean and ±SD of triplicate determinations are shown. (£) De novo and targeted (asterisk) identification of motifs under RXR peaks 
from ChlP-seq data using HOMER. "Target %" refers to the ratio of the peaks having the given motif, and "Bg %" refers to the ratio of 
a random background as described in the Supplemental Material. (F) The genome- wide distribution of RXR peaks relative to the TSS of 
the closest regulated genes identified by RNA-seq [P < 0.05; n = 823). 
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RXR peaks overlap with both PU.l and P300 to a large 
extent and show increased acetylation (Fig. 2A,B). Global 
analyses show that ligand activation of RXR leads to P300 
recruitment (Fig. 2B,C). These changes could be validated 
by ChlP-RT-qPCR (Fig. 2C). Importantly, global P300 
binding and recruitment could be detected and confirmed 
at >400 sites using three biological replicates, and no 
decrease could be observed (Supplemental Fig. S2A). As 
far as active histone marks are concerned, H4 acetylation 
and H3K4 dimethylation are increasing, while H3K27 
acetylation remains largely unchanged upon ligand acti- 
vation on the RXR-bound genomic regions (Fig. 2D). 

Next, we sought to determine whether some of the 
potential partner ligands had a detectable effect on the 
RXR cistrome. Therefore, we determined the RXR cis- 



tromes in the presence of a synthetic LXR (GW3965) or 
PPAR7 (RSG) ligand as the activators of the two major 
suspected heterodimeric partners present in macrophages. 
We found that neither ligand had a statistically significant 
effect on RXR's genome-wide distribution and/or affinity, 
applying a <0.1 false discovery rate (FDR) threshold (Fig. 
2E). Finally, we wanted to evaluate the possibility of 
whether an endogenous RXR ligand is masking a ligand- 
induced (re)distribution of the receptor. We reasoned that 
eliminating the effect of a presumed endogenous ligand by 
an RXR antagonist (LG1208) would reveal the true unli- 
ganded RXR cistrome. However, as shown in Figure 2E, 
the RXR cistrome does not change upon antagonist treat- 
ment either. These findings were confirmed using RT- 
qPCR (Supplemental Fig. S2B). These data suggest that the 
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Figure 2. The cistromic interactions of RXR, PU.l, P300, and active histone marks. (A) The intersection of cistromes as assessed by 
number of overlapping peaks is represented as a Venn diagram of RXR and PU.l. (B) Heat map representation of RXR, P300, H3K4me2, 
H3K27ac, and H4ac occupancies in 3-kb windows around the summit of the RXR peaks in the presence or absence of LG268. Read 
distribution was determined by HOMER, clustering was done by Gene Cluster 3.0 using centered correlation similarity metric with single 
linkage clustering method, and heat maps were created by Java TreeView in log 2 scale. (C) Read distribution of P300 on RXR peaks in the 
presence or absence of LG268. P300 binding was confirmed on the indicated individual enhancers using ChlP-RT-qPCR. Macrophages were 
treated with 100 nM LG268 for 1 h. The mean and ± SD of three biological replicates are shown. Asterisk represents significant difference at 
P < 0.05; n = 3. (D) Average read distribution of the indicated active histone marks on RXR peaks in the presence or absence of LG268 
determined by HOMER. (E) RXR enrichments of the significantly changing peaks in the presence of the indicated ligands. Cistrome was 
determined in the presence or absence of LG268, and the log2-normalized read numbers of the significantly changing peaks were plotted by 
DiffBind. The number of changing peaks and the statistical stringency applied is indicated below each plot. 



RXR cistrome is not changing upon short-term ligand 
activation and that the RXR selective ligand only increases 
the receptor's enrichment at its predetermined binding 
sites, which typically bind PU.l and recruit P300. Impor- 
tantly, the dominant ligands of relevant permissive hetero- 
dimers such as PPAR7 or LXR do not appear to have an 
effect on either distribution or affinity. 

Determination of liganded RXR regulated nascent 
RNA production 

The merging of the ChlP-seq and RNA-seq data sets 
proved to be insufficient to reliably identify the network 



of primary regulated genes and their enhancers. This is 
due to the fact that the physical location of the bound 
transcription factor cannot be linked to the regulated 
gene. Almost 90% of the RXR peaks are located outside 
of the 10-kb region upstream of the TSS of the closest 
regulated gene (Fig. IF), which is in good correlation with 
others 7 observations (Thurman et al. 2012). 

We reasoned that in order to link a liganded RXR- 
occupied enhancer to the corresponding regulated 
gene, one needs a high-resolution method, which can 
provide dynamic temporal information about transcrip- 
tion. The recently recognized fact that nascent RNA 
production can be detected on both the regulated gene 
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and the enhancer offered a plausible solution (Wang 
etal. 2011). 

Therefore, we decided to detect the dynamics of na- 
scent RNA production using GRO-seq. As shown in 
Figure 3A, we carried out a time-course experiment to 
determine the sites and dynamics of nascent RNA pro- 
duction upon LG268 activation, a typical example of 
which is shown for Abcal induction in Figure 3B. Impor- 
tantly, we noticed the presence and induction of nascent 
RNA transcripts at sites of enhancers, eTranscripts, and 
also short divergent transcripts (Core et al. 2008). We 
collectively call these enhancer RNAs (eRNAs). Nota- 
bly, if one aligns the detected GRO-seq activities with 
cistromic data, it is easy to recognize that short tran- 
scripts overlap with transcription factor-bound regions — in 



our case, PU.l, RXR, and P300 (Fig. 3C ; Supplemen- 
tal Fig. S3 A). A comprehensive analysis of GRO-seq 
active regions, including gene transcript and eRNA 
calling, was carried out by a unique algorithm (Supple- 
mental Fig. S3B). We determined that -20% of the 
mouse genome is transcribed (with at least 0.006 RPKM 
[reads per kilobase per million mapped reads] expression) 
in this macrophage cell type and identified 10,586 
known genes also marked by H3K4me3 and, all together, 
25,560 transcripts, which include eRNAs and other 
noncoding RNAs as well (Supplemental Fig. S3C). Im- 
portantly, we found 51,657 GRO-seq divergent sites 
(>0.2 RPKM expressed) characterized by divergent tran- 
scription with a distance <300 bases between the 5' ends 
of the transcripts on the two strands. These sites do not 
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Figure 3. Determination of nascent RNA production upon activation of RXR. [A) The scheme of GRO-seq experiments. Cells were 
treated for the indicated time periods and collected, nuclei were isolated, and run-on sequencing was performed. [B) The detected 
induction of GRO-seq activity on the Abcal locus. Nascent RNA produced around the Abcal locus is shown. Strand-specific coverage 
is represented by red and blue. (C) Genome browser view of the merge of GRO-seq, and ChlP-seq activity/peaks on the Abcal locus. 
Overlaps of PU.l, RXR, and P300 binding with divergent sites (GRO-seq-positive) are indicated with arrows. [D] Comparison of mRNA 
production dynamics using nascent (GRO-seq) and steady-state (RNA-seq and RT-qPCR) RNA determinations. Macrophages were 
treated with LG268 for the indicated time period, and specific gene expression was detected with the indicated method. A 
representative set of experiments is shown. Expression was normalized to the 0 time point. 
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necessarily show elongation, just initiation irrespective 
of transcribed genes and other longer transcripts. One 
could reasonably assume that this set of genomic regions 
contains all active promoters and enhancers. Using this 
approach, we determined changing (induced/repressed) 
gene transcripts and eRNAs as well. If one plots the 
changing levels of nascent RNA production of already 
established direct target genes, the dynamics of the 
changes are indicative of an immediate induction (Fig. 
3D). This also allows easy identification of additional 
direct targets. 

These data further reinforced the notion that the 
dynamically changing nascent RNA landscape provides 
clues about direct transcriptional regulation specific for 
RXR liganding. Furthermore, this allowed us to generate 
a list of genes most likely being directly regulated by 
liganded RXR (Supplemental Table SI). 

As far as the active enhancers are concerned, we 
classified the identified —5200 RXR-binding sites into 
two categories using an algorithm (pipeline depicted in 
Supplemental Fig. S4A): (1) enhancers that overlap with 
a GRO-seq-positive region (divergent site,- 2781) and (2) 
enhancers that do not overlap with GRO-seq activity 
(2425), as shown in Figure 4, A and B). This classification 
showed an enrichment of RXR peaks within the proxim- 
ity of the closest TSS (<10 kb from the TSS), suggesting 
that our assumption was likely correct and that we were 
identifying the functional/active binding regions and 
efficiently separating those from the silent/nonactive/ 
parking ones. The distribution of GRO-seq-positive and 
GRO-seq-negative RXR-bound regions relative to the 
nearest TSS shows a bias toward upstream regions in cases 
of the positive ones. Most interestingly, the enriched 
motif distribution under these peaks also showed charac- 
teristic differences. The GRO-seq-positive RXR sites 
showed an enrichment for DR1 and DR4 (binding sites of 
RXR:PPAR and RXR:LXR heterodimers, respectively) as 
well as AP-1 sites when compared with the GRO-seq- 
negative sites (Fig. 4B). 

GRO-seq-positive RXR-binding sites are characterized 
by an increased H3K27ac, H4ac, and also H3K4me2 
histone marks when compared with negative ones (Fig. 
4C,D). These provided further support to the notion that 
these sites are indeed functionally distinct and likely 
represent active enhancers. 

A closer look at the up-regulated GRO-seq-positive 
regions revealed that those that overlap with RXR peaks 
overlap with a significant fraction of PU. 1 regions as well 
(Fig. 5 A, top). The regions of down-regulated GRO-seq 
activity have much fewer RXR peaks but show overlap 
with PU.l regions to a larger degree as well (Fig. 5 A, 
bottom). These latter analyses showed that there are very 
few liganded RXR-occupied negative binding sites or 
silencers, suggesting also that liganded RXR is predomi- 
nantly a transcriptional activator. Next, we matched up 
the divergent GRO-seq sites with RXR sites (Fig. 5B ; 
Supplemental Fig. S5A). The identified 51,657 divergent 
GRO-seq sites likely contain all TSSs, short divergent 
transcripts, and eRNAs. More than 3300 of these change 
upon RXR ligand activation, and 718 overlap with RXR- 



binding sites as well (Fig. 5A,C). We identified 252 
regulated genes to which we could assign 414 RXR-bound 
regions, 387 enhancers, and 27 silencers using the criteria 
set by us (Fig. 5B,C; Supplemental Fig. S4A). 

The distribution of the various motifs within these 
potentially active enhancers show that —50% can be 
classified as a nuclear receptor-binding site repeat, an- 
other 20% is PU. 1 only, and a quarter is other or unknown 
(Fig. 5D). The number of enhancers per gene ranges from 
one to nine (Fig. 5E). The distribution of these putative 
enhancers relative to the TSS is fairly symmetrical, and 
they span long distances, with only 16.54% found within 
10 kb upstream of the TSS (Fig. 5F). The identified 387 
enhancers have a pronounced increase in RXR occu- 
pancy, P300 recruitment, and H4ac, while H3K4me2 
and H3K27ac are not different upon ligand treatment 
(Supplemental Fig. S4B). 

Taking these together, we ( 1 ) uncovered the core active 
enhancer network operated by liganded RXR; (2) se- 
parated these sites from nonfunctional, silent/parking 
RXR-binding sites,- and (3) paired the active enhancer 
network with the regulated genes using a set of criteria. 

Functional validation of novel distant and long-range 
enhancers 

Before we embarked on molecular validation of the 
enhancers, we intended to assess the contribution of the 
activated partner receptors to the induction of various 
RXR-regulated genes and eRNAs also. Therefore, we de- 
termined the induction of steady-state RNA levels of 
selected genes. We found that Vegfa and Hbegf were 
induced primarily by RXR liganding with some activity 
by RAR, Tgm2 and Ccl6 were both induced by RXR and 
RAR liganding, Abcgl was induced by RXR and LXR, and 
Angptl4 was induced by the RXR and PPAR7 ligands 
(Supplemental Fig. S5). 

Interestingly, primarily RXR-induced eRNA produc- 
tion could be detected on an enhancer assigned to Vegfa 
or Tgm2, while an enhancer of Abcgl also showed robust 
LXR ligand activation, as expected (Fig. 5G). These data 
suggested that the eRNAs can be easily validated and 
show ligand induction similar to the regulated genes and 
therefore most likely are linked. 

For functional validation, we chose 45 newly identified 
regulatory regions assigned to 23 genes (Supplemental 
Table SI). We selected 30 distant enhancers (>10 kb from 
the TSS), which earlier studies performing traditional 
enhancer analyses (Schwartz et al. 2000) would likely 
have missed, including one of ours (Nagy et al. 1996). In 
addition, we sought to validate 15 more proximal en- 
hancers (<10 kb from the TSS). 

Next, we cloned the 45 putative cis -regulatory ele- 
ments (35 enhancers and 10 silencers) in front of a lucif- 
erase reporter gene and measured their transcriptional 
activity in COS 1 cells in the presence of combinations of 
receptors and ligands. We clustered the enhancers based 
on the induction patterns obtained in the transient trans- 
fection-based reporter system, which is detailed in the 
Supplemental Material (Supplemental Fig. S6). These 
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Figure 4. Classification of RXR-binding sites based on GRO-seq activity. [A, top left) Distribution of the (GRO-seq + ) divergent sites 
detected based on genome-wide localization relative to the closest expressed transcripts (>3 kb including 3' overhang) determined by 
GRO-seq. Distribution of all RXR sites (top right), the RXR sites overlapping with GRO-seq activity [bottom left), and the RXR sites not 
overlapping with GRO-seq activity {bottom right) relative to the expressed transcripts. [B) Distribution of GRO-seq-positive (blue) and 
GRO-seq-negative (red) RXR peaks compared with the TSSs of the closest expressed transcripts (defined as in A). Columns represent the 
peak number of the 1-kb distance bins. {Insert) De novo and targeted (black asterisk) identification of motifs under GRO-seq-positive 
(left side) and GRO-seq-negative [right side) RXR peaks determined using HOMER. Differentially represented motifs are marked by red 
asterisks. "Target %" refers to the ratio of the peaks having the given motif, and "Bg %" refers to the ratio of a random background as 
described in the Supplemental Material. (C) Heat map representation of RXR, P300, H3K4me2, H3K27ac, and H4ac occupancies in 3-kb 
windows around the summit of the GRO-seq-positive and GRO-seq-negative RXR peaks in the presence or absence of LG268. Read 
distribution was determined by HOMER, clustering was done by Gene Cluster 3.0 using centered correlation similarity metric with 
single linkage clustering method, and heat maps were created by Java TreeView in log 2 scale. [D) Read distribution of the indicated 
histone marks relative to the GRO-seq-positive and GRO-seq-negative RXR cistromes determined by HOMER. 



analyses show that the identified and tested enhancers 
are RXR-specific enhancers but also complex and very 
versatile, allowing combinatorial and context-dependent 
regulation of the genes by distinct dimers. The context is 



determined by the receptor expression profile and prob- 
ably also lineage-specific transcription factors and co- 
factors. In addition, a single gene is likely to have distinct 
types of enhancers. 
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Figure 5. Identification of RXR-regulated enhancers and linking those to regulated genes. [A] Overlap of sites of changing GRO-seq 
activity and cistromes of liganded RXR and PU.l: increasing GRO-seq activity [top] and decreasing GRO-seq activity (bottom). [B) 
Scheme of annotation criteria on the Abcal locus. This includes increased (or decreased) nascent transcript production of a gene and 
increased (or decreased) divergent transcription of a region occupied by RXR. (C) Flow of annotation of GRO-seq-positive RXR-binding 
regions with changing enhancer transcription to the closest gene with changing expression. See also Supplemental Figure S4A. 
(D) Motif distribution of the annotated 387 RXR-bound enhancers determined by HOMER and fuzznuc as described in the 
Supplemental Material. Note that additional PU.l sites can be found in conjunction with other motifs, which are not included here. 
(£) Distribution of the number of enhancers assigned to the closest gene; the Y-axis shows the number of genes having the number of 
enhancers (X-axis). [F] Distribution of enhancers relative to TSSs of the identified directly regulated genes. (G) RT-qPCR analysis of 
Vegfa, Abcgl, and Tgm2 eRNAs. Differentiated macrophages were treated with LG268, RSG, GW3965, and AM580 for 1 h, and eRNAs 
were determined. Data represent mean and ±SD, with expression normalized to RplpO. 
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The genome architectural context of RXR-regulated 
gene expression 

Once we determined the RXR-regulated genes and linked 
at least a subset of them to their enhancers (some of them 
are long-range ones), we wanted to get some insight into 
whether these regulatory units are confined to known 
features of genome architecture. 

Therefore, we decided to explore the genome archi- 
tectural context of RXR signaling, asking the question of 
whether RXR signaling is confined by known structural 
features of genome organization and/or contributes to 
reorganization of these. We determined the cistromes of 
CTCF and a member of the cohesin complex (RAD21 ) in 
the absence or presence of ligand-activated RXR. We 
found -30,000 peaks for CTCF and -24,500 for RAD21 
(Supplemental Fig. SIC); the latter showed some in- 
crease upon ligand activation. Next, we determined an 
overlap of —12,660 peaks representing similarly high 
CTCF and RAD21 occupancy (with peaks having 
a MACS2 score >15 and having less than threefold 
difference). We considered these as boundaries of func- 
tional domains, as suggested by others (Merkenschlager 
and Odom 2013; Sofueva et al. 2013). We identified 10,204 
such functional domains by pairing the CTCF and RAD21 
copeaks with the closest ones with similar scores (Supple- 
mental Material). 

The median length of such domains was 81.15 kb, 
while the average was 149.98 kb. The neighboring do- 
mains then were merged into active domains with 
<100-kb distances between them. We could designate 
almost 700 such domains, which have a median length of 
1.1 Mb and an average length of 1.88 Mb. If one overlays 
this architectural domain structure with active RXR- 
binding sites and regulated genes that we identified (Fig. 
5C), one can find that 203 out of the 252 (80%) identified 
RXR-regulated genes along with their enhancers fall on 
such functional domains. In addition, one can identify 
— 1 1 13 CTCF genomic binding regions on which RAD21 
(cohesin) binding is enhanced and 128 regions on which 
RAD21 binding is reduced upon RXR ligand treatment 
(Fig. 6 A, left). Moreover, RAD21 (cohesin) binding is more 
enhanced on active RXR enhancers (GRO-seq + ones) 
when compared with nonactive ones (GRO-seq - ones) 
(Fig. 6A, middle and right). Importantly, 40% (84 out of 
203) of regulated genes have induced RAD21 on their 
enhancer and/or CTCF-binding sites. Collectively, these 
data suggest that most if not all RXR enhancers act inside 
functional domains and that activation of the receptor 
contributes to formation of such domains by stabilizing 
the genomic architecture and, in some cases, even in- 
ducing the binding of enhancers to promoters. Next, we 
validated molecularly the interactions in such functional 
domains and the impact of ligand on these for long-range 
enhancers of Abcgl, Vegfa, and Tgm2 using RT-qPCR- 
based 3C (Fig. 6B). These data clearly documented that 
the enhancers identified by the combination of RNA-seq, 
ChlP-seq, and GRO-seq act in functional domains, loop to 
the promoter, and can be readily functionally validated 
using transient transfection and 3C as well. 



Finally, we asked the question of whether any of these 
enhancers communicate with other cis-regulatory ele- 
ments or functional domains in the genome. Thus, we 
carried out 3C-seq using pairs of baits located in or close 
to these regions. We could detect proximal and also long- 
range interactions (for details of the analyses, see the 
Supplemental Material). At the Vegfa locus, we could 
detect interactions between the distant enhancer and the 
neighboring enhancers as well as the intronic region with 
remarkable specificity (Fig. 6C). Similar results were 
obtained in the case of the Tgm2 and Abcgl loci (Supple- 
mental Fig. S7A,B). Importantly, we also detected in- 
terchromosomal interactions with much less frequency 
(at least 50-fold less) though. In order to provide statistical 
context to these findings, we compared the interaction 
frequency of a given bait with inactive topological do- 
mains, active domains devoid of RXR-regulated regions, 
and active domains with RXR-regulated regions. As shown 
in Figure 6D in the case of Abcgl and Vegfa, there is 
significant difference between the frequencies of such 
interactions. Similar results were obtained with additional 
enhancer of these and the Tgm2 gene. This suggests that 
active RXR enhancers interact with other active genomic 
regions and with an even higher likelihood with other 
RXR-regulated ones. These formally suggest that the 
active RXR enhancers form an interchromosomal hub or 
network. 

Determination of the impact of RXR activation 
on the angiogenic capacity of macrophages 

Finally, we wanted to see whether some of the identified 
novel transcriptional pathways could be validated biolog- 
ically. The functional annotation of the genes controlled 
by liganded RXR enhancers assigned a large number of 
genes into the angiogenesis category (Supplemental Fig. 
S4C). These include Vegfa, Hbegf, Litaf, and Hipk2. 
Therefore, we decided to functionally test this activity. 

First, we determined the excreted VEGFa protein levels 
from cell supernatants (Fig. 7A) using different ligands 
and could show that RXR induced this protein, as was 
suggested by our gene expression measurements as well 
(Supplemental Fig. S5). Using RT-qPCRs, we could show 
the RXR requirement using RXRa/p double-knockout 
macrophages (Fig. 7B). To test whether activation of RXR 
promoted angiogenesis in an in vivo relevant setting, 
we performed a chorioallantoic membrane (CAM) assay 
using macrophages pretreated with RXR agonist. These 
cells showed significantly increased angiogenic activity, 
which was not detectable in RXRa/p double-knockout 
macrophages (Fig. 7C). These data suggest that RXR ac- 
tivation can program macrophages toward a distinct cell- 
autonomous angiogenic phenotype (Fig. 7D,E). 

Discussion 

The ultimate goal of studying the function of a particular 
transcription factor is to discern its primary biological 
action in the entire genome comprehensively. In princi- 
ple, this should be done by defining the genomic binding 
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Figure 6. Functional characterization of short- and long-range liganded RXR-regulated enhancers. (A) Changing RAD21 peaks in the 
presence of 100 nM LG268 on CTCF-binding sites (left) and GRO-seq-positive {middle) and GRO-seq-negative {right) RXR-binding 
sites. The log 2 -normalized read numbers of the significantly changing peaks are plotted. The number of changing peaks and the 
statistical stringency applied are indicated below each plot. (B) 3C-RT-qPCR measurements in the presence or absence of LG268 on the 
Tgm2, Abcgl, and Vegfa loci. Constant primers were designed to the enhancers and are represented by black columns (anchor points). 
Red circles represent the enhancers, purple ones are the promoter regions, TSSs are depicted with green, and intronic regions are 
highlighted by cyan circles. The mean and ±SD of triplicate determinations are shown. Representative examples are shown for each 
locus. (C) Genome browser view of the Vegfa locus containing the proximal interacting regions of the intergenic (Bl) and intronic (B2) 
baits and the loop predictions generated based on CTCF/RAD21-cobound regions. 3C-seq heat map (green scale) was made as described 
in the Supplemental Material. Asterisks show the sites of the specific baits. GRO-seq and ChlP-seq for the indicated factors are shown. 
Black arrowheads and gray dashed lines indicate the predicted domain borders. (D) Box plot representation of the distribution of 
interaction frequency of Abcgl Bl enhancer (chrl7: 31,172,998-31,173,656) and Vegfa Bl enhancer (chrl7: 45,890,060-45,890,829) 
determined by 3C-seq. The interchromosomal interactions of these enhancers were determined as described in the Supplemental 
Material. The detected interactions were mapped onto 1-Mb fragments covering the mouse genome. Interaction frequency was 
determined by expressing the read number per 1 Mb normalized to 1000 reads. Genome regions (inactive, active, and RXR target) were 
qualified as described in the Supplemental Material. Unpaired two-tailed t-test analysis was used to determine significant differences. 
Asterisk represents significant difference at P < 0.0001. 
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Figure 7. RXR activation leads to increased angiogenic activity in macrophages. (A) VEGFa protein levels determined from 
supernatants of macrophages treated with LG268 for the indicated periods of time determined by ELISA. The mean and ±SD of 
triplicate determinations are shown. Asterisk represents significant difference at P < 0.05; n = 3. [B] Vegfa mRNA levels determined 
using RT-qPCR in wild- type (WT) and RXRa/p knockout (KO) macrophages. The mean and ±SD of duplicate determinations are 
shown. (C) Macrophage angiogenesis activity was determined using CAM assay in the presence or absence of LG268 using wild-type 
and RXRa/p knockout macrophages. Representative images {left) and bar chart of quantification {right) are shown. The mean and ±SD 
of three biological replicates are shown. Asterisk represents significant difference at P < 0.05; n = 3. {D) The molecular mechanism 
through which Vegfa is regulated in macrophages. Vegfa harbors a set of very-long-range enhancers. The enhancers are marked by the 
lineage-determining transcription factor PU. 1 and are able to recruit RXR presumably as a heterodimer with RAR or in a homodimeric 
form. Upon ligand stimulation, P300 is recruited, and eRNA production is highly increased. These elements can get into the close 
proximity of the gene promoter by looping, which is enhanced in the presence of the RXR activator (LG268) and leads to increased gene 
expression. This regulatory unit is bordered by CTCF/cohesin (RAD21) cobinding, which is thought to contribute to the topological 
domain structure of mammalian genomes and stabilizes chromatin loops. (£) Schematic representation of the various nuclear receptor 
dimers and their biological significance in macrophages. PPAR/RXR are known to inhibit the inflammatory response upon activation 
by either side as a permissive heterodimer, while the LXR/RXR heterodimer regulates inflammatory response, cholesterol metabolism, 
and triglyceride biosynthesis. Our results shed light on the effect of putative, permissive, RAR/RXR, or possibly RXR/RXR dimers in 
the angiogenesis program of macrophages, most probably through the regulation of Vegfa and/or Hbegf. 
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regions and determining the genes proximally regulated 
in a given cell type. Using computational tools, one can 
identify potential RXR-binding sites on the scale of several 
times 100,000 in the mouse genome. When all of the 
ChlP-seq data available from fat, from liver, and now in 
macrophages are combined, the combined number of 
genomic regions bound by RXR is a few times 10,000. 
Still, in macrophages, the number of genomic regions is 
—5200, but the number of regulated genes by liganding the 
receptor is a few times 100. This is further complicated by 
the fact that primary targets are hard to reliably identify by 
simply measuring steady-state mRNA levels. So the 
question remains: How can one find the regulated genes 
and the linked receptor-bound enhancers? We attempted 
to do this by integrating genome-wide analyses. 

Integration of genome-wide localization and nascent 
RNA production data identifies active enhancers 

A key premise of the work presented here is that genome- 
wide localization studies, predominantly ChlP-seq exper- 
iments, can be intersected with data from nascent RNA 
determinations such as GRO-seq time courses, which 
then enables one to comprehensively annotate cistromes 
and identify active enhancers (Hah et al. 2013; Lam et al. 
2013). This relies on the assumption that in vitro (nuclear 
run-on) determination of the activity of RNA polymer- 
ases (RNAPs) is a reliable indicator of transcriptional 
activation. This has been proven in other systems and 
leads to the identification of active enhancers in other 
cell types (Koch and Andrau 201 1; Wang et al. 201 1; Bonn 
et al. 2012), although systematic comparisons or cross- 
validation with other methods such as subcellular RNA 
fractionation (Bhatt et al. 2012) has not been done. The 
usage of GRO-seq for the determination of direct tran- 
scriptional responses has two major advantages: (1) The 
dynamics of the nascent RNA levels depends only on the 
rate of RNAP activity, and therefore it is matching 
the expected time course of a directly regulated gene. 
Our data are clearly demonstrating this because known 
and established direct target genes (Abcal, Abcgl, 
Angptl4, and Tgm2) as well as newly identified ones 
(i.e., Vegfa and Hbegf) show an immediate induction 
when assessed by GRO-seq rather than a complex, often 
delayed one determined by RNA-seq or RT-qPCR (Fig. 
3D). (2) Enhancers contain engaged RNAPs, and their 
activity can be revealed by GRO-seq in the form of 
a typically few-hundred-base-pair-long divergent tran- 
script, termed eRNA. The existence of such transcripts 
has been shown before (Kim et al. 2010; Wang et al. 201 1). 
Although the mechanism of their production and role in 
transcription is far from being clarified, it can be reason- 
ably assumed that they represent active enhancers, as 
the presence of an engaged RNAP strongly indicates it. 
Furthermore, if the changes of nascent RNA production 
(increase or decrease) at these sites show correlation with 
the stimulating signal, it is likely to be associated with or 
the consequence of the activity of the signal-specific 
transcription factor. A key limitation of today's technology 
is the lack of high-throughput validation of the identified 



enhancers. Nonetheless, our high rate of success in vali- 
dating such enhancers with low-throughput enhancer trap 
approaches, 3C-seq and 3C-RT-qPCR, suggests that signal- 
specific changes in eRNA production can be used to filter 
active enhancers from the general pool of genomic binding 
regions. Altogether, this also means that this combined 
approach can serve as proof of concept and a model to 
tackle similar problems with other signal-specific tran- 
scription factors. A remarkable feature of the approach is 
that it allowed the identification of active enhancers in 
spite of the fact that these represent a fraction of the 
binding regions (not more than 10%-15% by using our 
stringent criteria), including very-long-range enhancers. It 
is also interesting that the enhancers 7 distribution relative 
to the TSS appears to be rather symmetrical instead of 
being much biased to the upstream regions, suggesting 
that one needs to look in both directions and also far away 
to identify the true enhancer controlling a particular gene 
in a given cellular context. A cautionary conclusion of this 
study is that it is very likely that many of the enhancers 
identified with less comprehensive methods might need to 
be revisited. 

RXR colocalizes with PU.l and recruits P300 
as a cof actor to its binding sites 

Our data showed that RXR has >5200 binding regions in 
macrophages. These peaks are likely to be RXRa -bound 
because this is the dominant subtype present in macro- 
phages. However the used antibody is pan-RXR-specific 
and would recognize all three subtypes. Most of these 
sites contain nuclear receptor half sites (AGGTCA), 
suggesting that the localization of RXR uses nuclear 
receptor-mediated direct DNA binding. In addition, two 
independent lines of evidence suggest that there is an 
intimate relationship between the lineage-determining 
factor PU.l and RXR: (1) 45% of the RXR peaks contain 
PU.l sites, and (2) there is a significant overlap between 
RXR and PU. 1 cistromes, as more than two out of three of 
RXR-binding regions overlap with PU.l peaks. This is in 
agreement with our anticipation based on previous re- 
ports by Lazar and colleagues (Lefterova et al. 2010) and 
Liu and colleagues (Pott et al. 2012), who suggested that 
binding regions for PPAR7, a heterodimerization partner 
of RXR, colocalizes with PU. 1 . This finding supports the 
pioneering or bookmarking factor concept put forward 
by several investigators (for a review see, Zentner and 
Scacheri 2012) to explain various interactions and ge- 
nome-wide colocalization between lineage-specific and 
signal-specific transcription factors, suggesting that the 
lineage-specific factor opens up or marks particular re- 
gions in the genome, which then allows or even facilitates 
the binding of the signal-specific factors. Although the 
sequence of events or the mechanisms of such interac- 
tions are unknown, our data are fully compatible with 
such a scenario. An additional line of evidence support- 
ing this is that there is only a 14% overlap between our 
RXR-binding regions and the ones found in 3T3L1 cells 
differentiating into adipocytes (Nielsen et al. 2008). Our 
work further extended these studies and suggests that 
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many if not all RXR heterodimers at least in part colo- 
calize with PU.l in macrophages and that upon ligand 
treatment, PU. 1 is at least partially released and therefore 
might not be required for further binding or transcrip- 
tional activity. 

Another unexpected aspect of the cistromic interac- 
tions identified is that P300 binding is dynamically 
recruited upon the RXR ligand treatment to most of the 
active RXR sites. Two conclusions can be drawn from 
these findings: (1) P300 is major cofactor of mediating 
transcriptional activity by RXR in BMDMs in the steady 
state. (2) Liganded RXR-regulated transcription complexes 
use P300 as one and potentially the major cofactor, and 
thus P300 is likely to be responsible for the acetylation of 
histone H4K5/8. This is consistent with the initially 
proposed role for this protein in nuclear receptor signaling 
(Chakravarti et al. 1996) and with the more recent findings 
(Jin et al. 201 1) but also suggests that P300 is likely to be 
specific for this signaling pathway in this cell type. 

A key issue with signal-specific transcription factors is 
whether the signal contributes to the (re (distribution 
of the particular transcription factor. In case of the 
RXR-containing receptor dimers, the accepted view is 
that RXRs bind to the genome in both the absence and 
presence of ligand (Boergesen et al. 2012). Our findings 
further support this view now by adding a genome-wide 
perspective and show that the only effect RXR ligand has 
on the genome-wide distribution of the receptor is that 
the enrichment on the preformed binding sites increases. 
These data obviously have to be interpreted within the 
context of ChIP experiments, meaning that the sum of 
all binding in all cells is detected in the time frame of 
the 40-min cross-linking used. This method does not 
allow the construction of a more dynamic and/or higher- 
resolution picture of the receptors 7 activation and mobil- 
ity. Methods using shorter time resolution show a much 
more dynamic behavior of the receptor though (Brazda 
et al. 2014). Therefore, our interpretation of the data is 
that RXRs 7 genome-wide localization is determined by its 
own DNA-binding capacity, which allows it to find 
preformed sites in the genome to which it can bind, and 
this is facilitated by additional factors such as a lineage- 
specific factor,- i.e., PU.l. The increased affinity that is 
observed is either the reflection of more cells being 
involved in the response or a higher affinity for the 
binding site. However, the experimental approaches used 
are not able to determine the relative contribution of 
these two mechanisms. We also excluded the possibility 
that an endogenous ligand plays a major role in directing 
RXRs to its genomic binding sites by using an antagonist. 
This behavior is in stark contrast to that of steroid hor- 
mone receptors such as GR or estrogen receptor (ER), 
whose genomic binding is dictated by the addition of ligand 
(Carroll et al. 2005; Welboren et al. 2009). As far as the 
genome architectural context of the RXR-regulated en- 
hancers and the regulated genes are concerned, we could 
show that the vast majority of enhancers and their regu- 
lated genes are confined to functional domains bordered by 
CTCF/cohesin (RAD21) complexes, suggesting that their 
activity localized to these loops (Supplemental Fig. S8). 



In addition, the rearrangement of such loops is moderately 
though but impacted by ligand activation. Curiously 
though, our 3C-seq analyses also revealed that the func- 
tional domains and the enhancers within them interact 
with each other even on different chromosomes. The func- 
tional significance of such interactions is not clear but 
might suggest the existence of RXR- or partner-specific 
transcription foci or factories inside the nucleus. 

Characteristics of the liganded RXR-operated 
enhancer network 

Our data show, as expected, that the two major permis- 
sive heterodimers in BMDMs are RXR:LXR and, to 
a lesser extent, RXR:PPAR, as shown by the ready 
induction of their established target genes, Angptl4, 
Abcal, and Abcgl, and the motifs identified within 
their RXR peaks being RXR:PPAR-specific (DR1) and 
RXR:LXR-specific ones (DR4). However, our analyses 
also uncovered a set of regulated genes, including Tgm2, 
and novel ones, such as Vegfa and Hbegf, which could 
not be induced efficiently by ligands for permissive 
heterodimers, only the RXR-selective LG268 or ligands 
activating the RAR receptor. This raises the possibility 
of the existence of permissive RXR:RAR heterodimers 
or a complex regulatory mechanism/enhancer allowing 
activation by either RAR or RXR-activating ligands 
(Fig. 7D). 

Underpinning this complexity, the enhancers identi- 
fied for the regulated genes show a large degree of func- 
tional versatility and can be grouped into four broad 
categories. Cluster A requires the presence of liganded 
RXR; another (cluster B) is best activated if both RXR 
and one of its partners is expressed; a third (cluster C) is 
mediating RXR-specific signaling, provided RXR is 
expressed at a significantly higher level than its partners,- 
and the enhancers in the fourth (cluster D) work with 
multiple combinations of dimers (Supplemental Fig. S4). 
A further level of complexity is that a given gene (i.e., 
Vegfa or Abcal) has enhancers from more than one 
category. 

All of these support that RXR has a unique genomic 
effect that cannot be recapitulated by any ligand or 
combination of ligands. This is mechanistically served 
by activation of permissive heterodimers such as LXR:RXR 
and PPAR:RXR and activation of enhancers inducible 
by RAR as well. In addition, some of the eTranscript 
inductions appear to be specific for RXR. 

This and the facts that, using our most stringent 
criteria, we identified only 226 induced genes, and only 
7.4% of the identified binding sites are active enhancers 
suggest that activation of RXR is tightly controlled in 
this cell type and leads to distinct and selective and not 
pleiotropic gene expression. Unexpectedly, a network of 
enhancers that could be linked to a set of genes — includ- 
ing Vegfa, Hbegf, Hipk2, Litaf, Cxcl2, and Foxo3, which 
have been clearly linked to angiogenesis — was also 
revealed. 

The regulation of Vegfa by RXR is intriguing because it 
appears to use a set of very-long-range enhancers (270 kb 
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downstream) and shows RXR/RAR specificity, including 
RXR ligand-induced looping (Fig. 7D). We presented here 
a novel integrated approach to identify functional en- 
hancers and link those to primarily regulated genes for 
a signal-dependent transcriptional factor, liganded RXR. 
This integrated approach revealed that ligand stimulation 
of RXR activates only a small fraction of the DNA-bound 
molecules confined by CTCF/cohesin-delimited func- 
tional domains (Supplemental Fig. S8) and leads to a dis- 
tinct gene expression program, which results in increased 
angiogenic potential and might be a valid macrophage 
reprograming/therapeutic target (Fig. 7E). 

Materials and methods 

Materials 

The following ligands were used: LG268 and LG1208 were gifts 
from M. Leibowitz (Ligand Pharmaceuticals), RSG and AM580 
were from Sigma, and GW3965 was a gift from T.M. Wilson 
(GlaxoSmithKline). 

Differentiation of BMDMs 

Isolation and differentiation were completed as described earlier 
(Barish et al. 2005). 

RNA-seq 

The RNA-seq library was prepared from two biological replicates 
by using a TruSeq RNA sample preparation kit (Illumina) 
according to the manufacturer's protocol. Analysis was carried 
out as described in the Supplemental Material. 

ChIP 

ChIP was performed as previously described (Barish et al. 2010), 
with minor modifications. Libraries were prepared by Ovation 
Ultralow Library Systems (Nugen) from two biological replicates 
according to the manufacturer's instructions. ChlP-seq analysis 
and domain prediction were carried out as described in the 
Supplemental Material. 

GRO-seq 

GRO-seq and library preparation were performed as described 
earlier (Core et al. 2008; Hah et al. 2011), with limited modifi- 
cations. Libraries were generated from two biological replicates 
of BMDMs treated with 100 nM LG268. Analysis was carried out 
as described in the Supplemental Material. 

RT-qPCR 

RNA was isolated with Trizol reagent (Molecular Research 
Center). RNA was reverse-transcribed with Tetro reverse tran- 
scriptase (Bioline). Transcript quantification was performed by 
RT-qPCR reaction using SYBR Green dye (Diagenode). Tran- 
script levels were normalized to Ppia or RplpO. 

3C 

3C experiments were completed as described previously with 
minor modifications (Miele et al. 2006). The detailed protocol is 
available in the Supplemental Material. 



3C-seq 

Experiments were carried out as previously described (Stadhouders 
et al. 2013). For details, see the Supplemental Material. 

Transient trans fection 

Enhancer sequences were PCR-amplified from a BAC or genomic 
DNA and cloned into pUC18 HSV TK-LUC. Transient trans- 
fections were carried out as previously described (Szanto et al. 
2010). 

ELISA 

ELISA experiments were carried out according to the manufac- 
turer's instructions (R&D Systems). 

CAM assay 

CAM assays were performed as described previously (Movahedi 
et al. 2008). The detailed protocol is presented in the Supple- 
mental Material. 

Data access 

Sequencing data were submitted to Sequence Read Archive 
under accession number SRP0 19970. 
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